{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "REJqvG0z6Bkj"
   },
   "source": [
    "# Randomness\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Y1bnC_p-66WG"
   },
   "source": [
    "What type of problems can we solve with the help of random numbers?\n",
    "\n",
    "We can compute (potentially) complicated averages:\n",
    "\n",
    "*   How much my stock/option portfolio is going to be worth?\n",
    "*   What are my odds to win a certain competition?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ONUZol4M8J7y"
   },
   "source": [
    "### Random Number Generators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "A0z8frFY8S5H"
   },
   "source": [
    "*   Computers are deterministic - operations are reproducible\n",
    "*   How do we get random numbers out of a deterministic machine?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "ETXJlNHY8pIa"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the library `numpy.random` to generate random numbes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 51
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "ZLgxPkJD8wYZ",
    "outputId": "7c56d135-4201-4023-d99e-6fb6e93b9a91",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "np.random.rand(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to generate a random integer number over a given range, you can use\n",
    "\n",
    "`np.random.randint(low,high)`\n",
    "\n",
    "that returns a random integer from low (inclusive) to high (exclusive)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.random.randint(1,10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that if you use the library `random` to accomplish the same thing:\n",
    "\n",
    "`random.randint(low,high)`\n",
    "\n",
    "the function returns a random integer from low (inclusive) to high (**inclusive**)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.random.randint(1,10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generating many random numbers at one, using a numpy array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 357
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "W_INXEwi8kBU",
    "outputId": "c31cb2c2-fc8e-4439-fb80-923cc56b3f0c"
   },
   "outputs": [],
   "source": [
    "for x in range(0, 20):\n",
    "    numbers = np.random.rand(6)\n",
    "    #print(numbers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "fFZS3WxE871v"
   },
   "source": [
    "They all seem random correct? Let's try to fix something called **seed** using  \n",
    "np.random.seed(10)\n",
    "\n",
    "What do you observe? \n",
    "\n",
    "Let's see what this seed is..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "IuiZnPxm9dQc"
   },
   "source": [
    "## Pseudo-random Numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "WRW753Be9ibZ"
   },
   "source": [
    "* Numbers and sequences appear random, but they are actually reproducible\n",
    "* Great for algorithm developing and debugging\n",
    "* How truly \"random\" are these numbers?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "M828HMrX94T9"
   },
   "source": [
    "## Linear congruential generator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Y7Eva39j-rkS"
   },
   "source": [
    "Given the parameters $a$, $c$, $m$ and $s$, where $s$ is the seed value, this algorithm will generate a sequence of pseudo-random numbers:\n",
    "\n",
    "$x_o = s $\n",
    "\n",
    "$x_{n+1} = (a x_n + c) mod(m)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "gZYMuH2t_ZZ5"
   },
   "outputs": [],
   "source": [
    "s = 3 # seed\n",
    "a = 37 #10 # multiplier\n",
    "c = 2 # increment\n",
    "m = 19 # modulus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "ePCLntyl_kqN"
   },
   "outputs": [],
   "source": [
    "n = 60\n",
    "x = np.zeros(n)\n",
    "x[0] = s\n",
    "for i in range(1,n):\n",
    "    x[i] = (a * x[i-1] + c) % m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 282
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "jzE_vTlw_4te",
    "outputId": "4cf92953-d8db-4e10-d250-15e7edf66469"
   },
   "outputs": [],
   "source": [
    "plt.plot(x,'.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "P-iFLPo5AVTr"
   },
   "source": [
    "Notice there is a period, when numbers eventually start repeating.  One of the advantages of the LCG is that by using appropriate choice for the parameters, we can obtain known and long periods.\n",
    "\n",
    "Check here https://en.wikipedia.org/wiki/Linear_congruential_generator for a list of commonly used parameters of LCGs.\n",
    "\n",
    "Let's try using the parameters from \n",
    "[Numerical recipes](https://en.wikipedia.org/wiki/Numerical_Recipes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 293
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "6GF31MZ6-oQP",
    "outputId": "d14ece91-7508-43db-8a98-d25ce4a345db"
   },
   "outputs": [],
   "source": [
    "s = 8\n",
    "a = 1664525\n",
    "c = 1013904223\n",
    "m = 2**32\n",
    "\n",
    "n = 300\n",
    "x = np.zeros(n)\n",
    "x[0] = s\n",
    "for i in range(1,n):\n",
    "    x[i] = (a * x[i-1] + c) % m\n",
    "\n",
    "plt.plot(x,'.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JV4q29UYBv85"
   },
   "source": [
    "\"Good\" random number generators are efficient, have long peiods and are portable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "y4CBQZcbB-8w"
   },
   "source": [
    "# Random Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "gWfjziCcCffv"
   },
   "source": [
    "Think of a random variable $X$ as a function that maps the outcome of an unpredictable (random) processses to numerical quantities.\n",
    "\n",
    "For example:\n",
    "\n",
    "* $X$ = the face of a bread when it falls on the ground. The random value can abe the \"buttered\" side or the \"not buttered\" side\n",
    "* $X$ = value that appears on top of dice after each roll\n",
    "\n",
    "We don't have an exact number to represent these random processes, but we can get something that represents the **average** case. To do that, we need to know the likelihood of each individual value of $X$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "4t4D9HC9FRpR"
   },
   "source": [
    "### Coin toss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Ar-qGRL8FMDK"
   },
   "source": [
    "Random variable $X$: result of a coin toss\n",
    "\n",
    "In each toss, the random variable can take the values $x_1 = 0$ (tail) and $x_2 = 1$ (head), and each $x_i$ has probability $p_i = 0.5$. \n",
    "\n",
    "The **expected value** of a discrete random variable is defined as:\n",
    "\n",
    "$$ E(x) = \\sum_{i=1}^{m} p_i x_i $$\n",
    "\n",
    "Hence for a coin toss we have:\n",
    "\n",
    "$$ E(x) = 1(0.5) + 0(0.5) = 0.5 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RAy-VqBRPlLN"
   },
   "source": [
    "### Roll Dice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JEdBUyl9PlLS"
   },
   "source": [
    "Random variable $X$: value that appears on top of the dice after each roll\n",
    "\n",
    "In each toss, the random variable can take the values $x_i = 1,2,3,...,6$ and each $x_i$ has probability $p_i = 1/6$. \n",
    "\n",
    "The **expected value** of the discrete random variable is defined as:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "4qPwU0PMRGVW",
    "outputId": "0cd295ac-3aa5-4f0d-fc8d-1048cf1ff9c3"
   },
   "outputs": [],
   "source": [
    "E = 0\n",
    "for i in range(6):\n",
    "    E += (i+1)*(1/6)\n",
    "E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "SoAgVnqV57hs"
   },
   "source": [
    "#  Monte Carlo Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "eLrLy0W9otyv"
   },
   "source": [
    "Monte Carlo methods are algorithms that rely on repeated random sampling to approximate a desired quantity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "gZDK6ifupA0s"
   },
   "source": [
    "### Example 1) Simulating a coin toss experiment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "T0arakhyR_rE"
   },
   "source": [
    "We want to find the probability of heads when tossing a coin. We know the expected value is 0.5. Using Monte Carlo with N samples (here tosses), our estimate of the expected value is:\n",
    "\n",
    "$$E = \\frac{1}{N}\\sum_{i=1}^N f(x_i) = \\frac{1}{N}\\sum_{i=1}^N x_i$$ \n",
    "\n",
    "where $x_i = 1$ if the toss gives head."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "KpTgSEa2pKLa"
   },
   "source": [
    "Let's toss a \"fair\" coin N times and record the results for each toss.\n",
    "\n",
    "But first, how can we simulate one toss?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "odxGZRnXSmoM",
    "outputId": "37a00f28-37f9-4266-d4ad-9a3ddc6a396a"
   },
   "outputs": [],
   "source": [
    "toss = np.random.choice([0,1])\n",
    "print(toss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "U9ZYKgj0SSVW",
    "outputId": "7e95e3c0-5734-452c-8f90-e2fe5e840593"
   },
   "outputs": [],
   "source": [
    "N = 30 # number of samples (tosses)\n",
    "\n",
    "toss_list = []\n",
    "for i in range(N):\n",
    "    toss = np.random.choice([0,1])\n",
    "    toss_list.append(toss)\n",
    "\n",
    "np.array(toss_list).sum()/N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Efi5ZoVnrhZr"
   },
   "source": [
    "Note that if we run the code snippet above again, it is likely we will get a different result. What if we run this many times? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 333
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "bFTFeML5FGa0",
    "outputId": "be443173-472d-4c55-e217-4dd879d88ee7"
   },
   "outputs": [],
   "source": [
    "#clear\n",
    "N = 1000 # number of tosses\n",
    "M = 1000 # number of numerical experiments\n",
    "nheads = []\n",
    "for j in range(M):\n",
    "    toss_list = []\n",
    "    for i in range(N):\n",
    "        toss_list.append(np.random.choice([0,1]))\n",
    "    nheads.append( np.array(toss_list).sum()/N )\n",
    "nheads = np.array(nheads)\n",
    "\n",
    "plt.hist(nheads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.hist(nheads, bins=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "NYuuPvtXCDa-",
    "outputId": "fa9e8002-f59c-4c58-d35a-f31057b6b508"
   },
   "outputs": [],
   "source": [
    "print(nheads.mean(),nheads.std())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rM1o1S8HnT62"
   },
   "source": [
    "What happens when we increase the number of numerical experiments?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "9NZ-2H5LteBa"
   },
   "source": [
    "### Monte Carlo to approximate integrals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "QWpwUkiVvxq3"
   },
   "source": [
    "One of the most important applications of Monte Carlo methods is in estimating volumes and areas that are difficult to compute analytically. Without loss of generality we will first present Monte Carlo to approximate two-dimensional integrals. Nonetheless, Monte Carlo is a great method to solve high-dimensional problems. \n",
    "\n",
    "To approximate an integration\n",
    "\n",
    "$$ A = \\int_{x_1}^{x_2} \\int_{y_1}^{y_2} f(x,y) dx dy $$\n",
    "\n",
    "we sample points uniformily inside a domain $D = [x_1,x_2] \\times [y_1,y2]$, i.e. we let $X$ be a uniformily distributed random variable on $D$. \n",
    "\n",
    "\n",
    "Using Monte Carlo with N sample points, our estimate for the expected value (that a sample point is inside the circle) is:\n",
    "\n",
    "$$ S_N = \\frac{1}{N} \\sum_{i=1}^{N} f(X_i) $$\n",
    "\n",
    "which gives the approximate for the integral:\n",
    "\n",
    "$$ A_N = (x_2 - x_1)(y_2-y_1) \\frac{1}{N} \\sum_{i=1}^{N} f(X_i) $$\n",
    "\n",
    "Law of large numbers:\n",
    "\n",
    "as $N \\rightarrow \\infty$, the sample average $S_N$ converges the the expected value $E(X)$ and hence $A_N \\rightarrow A$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2) Approximate the area of a circle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Ycw4At-YtnHb"
   },
   "source": [
    "We will use Monte Carlo Method to approximate the area of a circle of radius R = 1.\n",
    "\n",
    "Let's start with a uniform distribution on the unit square  [0,1]×[0,1] . Create a 2D array samples of shape (2, N):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "VOcrb1ItuANJ"
   },
   "outputs": [],
   "source": [
    "N = 10**2\n",
    "samples = np.random.rand(2, N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xrgFHRgHuLXv"
   },
   "source": [
    "Scale the sample points \"samples\", so that we have a uniform distribution inside a square $[-1,1]\\times [-1,1]$. Calculate the distance from each sample point to the origin $(0,0)$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "CHci1eonuPtr"
   },
   "outputs": [],
   "source": [
    "xy = samples * 2 - 1.0 # scale sample points\n",
    "r = np.sqrt(xy[0, :]**2 + xy[1, :]**2)  # calculate radius"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(xy[0,:], xy[1,:], 'k.')\n",
    "plt.axis('equal')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "vb-GIrEAuisk"
   },
   "source": [
    "We then count how many of these points are inside the circle centered at the origin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "-IgDXBPDuWvL",
    "outputId": "5f6c7250-22f2-47a8-c14a-c6ef688f848d"
   },
   "outputs": [],
   "source": [
    "incircle = (r <= 1)\n",
    "count_incircle = incircle.sum()\n",
    "print(count_incircle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "KB2AZ1Htvnca"
   },
   "source": [
    "And the approximated value for the area is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "ifu3O7hKvmcB",
    "outputId": "0cb43347-740f-4de5-fc77-1547a44b96fa"
   },
   "outputs": [],
   "source": [
    "A_approx = (2*2) * (count_incircle)/N\n",
    "A_approx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "akngnIAlvLxm"
   },
   "source": [
    "We can assign different colors to the points inside the circle and plot (just for vizualization purposes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 279
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "KVbBFZz2u8ur",
    "outputId": "969d095a-e802-4f2d-aeef-5bba28537c65"
   },
   "outputs": [],
   "source": [
    "plt.plot(xy[0,np.where(incircle)[0]], xy[1,np.where(incircle)[0]], 'b.')\n",
    "plt.plot(xy[0,np.where(incircle==False)[0]], xy[1,np.where(incircle==False)[0]], 'r.')\n",
    "plt.axis('equal')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "W1Hcg1sLzuby"
   },
   "source": [
    "Combine all the relevant code above, so we can easily run this numerical experiment for different sample size N."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "1adLD05ivgXo",
    "outputId": "cc11d838-499f-4400-ae09-235763a87940"
   },
   "outputs": [],
   "source": [
    "#clear\n",
    "N = 10**2\n",
    "samples = np.random.rand(2, N)\n",
    "xy = samples * 2 - 1.0 # scale sample points\n",
    "r = np.sqrt(xy[0, :]**2 + xy[1, :]**2)  # calculate radius\n",
    "incircle = (r <= 1)\n",
    "count_incircle = incircle.sum()\n",
    "A_approx = (2*2) * (count_incircle)/N\n",
    "print(A_approx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "lZHO_7iY0QUM"
   },
   "source": [
    "Perform the same above, but now store the approximated area for different N, and plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "WW315PBL0DVU"
   },
   "outputs": [],
   "source": [
    "#clear\n",
    "N = 10**6\n",
    "samples = np.random.rand(2, N)\n",
    "xy = samples * 2 - 1.0 # scale sample points\n",
    "r = np.sqrt(xy[0, :]**2 + xy[1, :]**2)  # calculate radius\n",
    "incircle = (r <= 1)\n",
    "N_samples = np.arange(1,N+1)\n",
    "A_approx = 4 * incircle.cumsum() / N_samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JjPg1zwp0_D8"
   },
   "source": [
    "The approximated area is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "I0COCiQz1DAx",
    "outputId": "419551a5-a5ad-484b-cb10-2ace2d80363b"
   },
   "outputs": [],
   "source": [
    "#clear\n",
    "A_approx[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 282
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "f5P8RBO82H-b",
    "outputId": "28a8d28a-7128-4d65-da33-802aae137c84"
   },
   "outputs": [],
   "source": [
    "plt.plot(A_approx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qxdOIwO41CYj"
   },
   "source": [
    "Which as expected gives an approximation for the number $\\pi$, since the circle has radius 1. Let's plot the error of our approximation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 300
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "DBsy2SHj0v9C",
    "outputId": "652d44ab-5f18-433a-d48a-f7c270e944d8"
   },
   "outputs": [],
   "source": [
    "plt.loglog(N_samples, np.abs(A_approx - np.pi), '.')\n",
    "plt.xlabel('n')\n",
    "plt.ylabel('error')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 300
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "buaZ18zO1Yuh",
    "outputId": "ec1bd9d7-f1f1-4936-f2e4-d007b01a1dd8"
   },
   "outputs": [],
   "source": [
    "\n",
    "plt.loglog(N_samples, np.abs(A_approx - np.pi), '.')\n",
    "plt.xlabel('n')\n",
    "plt.ylabel('error')\n",
    "\n",
    "plt.loglog(N_samples, 1/N_samples**2, '.')\n",
    "plt.loglog(N_samples, 1/N_samples, 'r')\n",
    "plt.loglog(N_samples, 1/np.sqrt(N_samples), 'm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 333
    },
    "colab_type": "code",
    "collapsed": true,
    "id": "RSlwPXwE25pp",
    "outputId": "80f3269e-bb25-42e1-cee4-89aea08c0e12"
   },
   "outputs": [],
   "source": [
    "N = 10**3\n",
    "M = 1000\n",
    "\n",
    "A_list = []\n",
    "\n",
    "for i in range(M):\n",
    "    samples = np.random.rand(2, N)\n",
    "    xy = samples * 2 - 1.0 # scale sample points\n",
    "    r = np.sqrt(xy[0, :]**2 + xy[1, :]**2)  # calculate radius\n",
    "    incircle = (r <= 1)\n",
    "    count_incircle = incircle.sum()\n",
    "    A_list.append( (2*2) * (count_incircle)/N )\n",
    "\n",
    "A_array = np.array(A_list)\n",
    "\n",
    "plt.hist(A_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nMMNj5E-5Ans"
   },
   "source": [
    "### Example 3) Approximating probabilities of a Poker game"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "p77MexiL5V0o"
   },
   "source": [
    "<img src=\"figures/holdem.png\" alt=\"drawing\" width=\"400\"/>\n",
    "\n",
    "What is the probability of winning a hand of Texas Holdem for a given starting hand? This is a non-trivial problem to solve analytically, specially if you include many players and the chances that a player will fold (give-up) before the end of a round of cards. Let's use Monte Carlo Methods to approximate these probabilities!\n",
    "\n",
    "Assumptions:\n",
    "- You are playing against only one player (the opponent)\n",
    "- Both players stay until the end of a game (all 5 dealer cards), so players do not fold.\n",
    "\n",
    "Monte Carlo simulation: for a given \"starting hand\" that we wish to obtain the winning probability, generate a set of N games and use the Texas Holdem rules to decide who wins (starting hand, opponent or there is a tie).\n",
    "\n",
    "A game is a set of 7 random cards:  2 cards for opponent + 5 cards for the dealer.\n",
    "\n",
    "For example, suppose the starting hand is 5 of clubs and 4 of diamonds. You perform N=1000 games, and counts 350 wins, 590 losses and 60 ties. Your numerical experiment estimated a probability of winning of 35%.\n",
    "\n",
    "**This is your first MP!**\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 4) Calculating a Volume of Intersection\n",
    "\n",
    "In this exercise, we will use Monte Carlo integration to compute a volume of intersection between two cylinders. This integral is possible to compute analytically, so we can compare our answer to the true result and see how accurate it is.\n",
    "\n",
    "The solid common to two right circular cylinders of equal radii intersecting at right angles is called the Steinmetz solid.\n",
    "\n",
    "Two cylinders intersecting at right angles are called a bicylinder or mouhefanggai (Chinese for \"two square umbrellas\").\n",
    "\n",
    "![](figures/steinmetz.JPG)\n",
    "\n",
    "http://mathworld.wolfram.com/SteinmetzSolid.html\n",
    "\n",
    "To help you check if you are going in the right direction, you can copy the functions you define here inside PrairieLearn.\n",
    "\n",
    "https://prairielearn.engr.illinois.edu/pl/course_instance/52088/assessments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### a) Write a function that will determine if a given point is inside both cylinders\n",
    "\n",
    "Write the function `insideCylinders` that given a NumPy array representing some arbitrary point in a 3-dimensional space returns `true` if the point is inside both cylinders. Assume the solid is centered at the origin, the cylinders are along the $x$ and $z$ axes and both have radius $r$\n",
    "\n",
    "```python\n",
    "def insideCylinders(pos,r):\n",
    "    # pos = np.array([x,y,z])\n",
    "    # r = radius of the cylinders\n",
    "    return bool\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### b) Write a function to evaluate the probability the point is inside the given volume\n",
    "The function `prob_inside_volume` should take as argument the number of random points N.\n",
    "\n",
    "The function generate N random points inside a box around the intersection of the cylinders, and uses the function `insideCylinders` to determine if the point is inside the cylinders or not. Recall that these random points should be generated in a form of a NumPy array.\n",
    "\n",
    "Track the number of points $C$ that fall inside both cylinders. Return the ratio $C/N$ as a floating point number.\n",
    "\n",
    "```python\n",
    "def prob_inside_volume(N,r):\n",
    "    # N = number of sample points\n",
    "    # r = radius of the cylinders\n",
    "    return float\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### c) Use the ratio $\\frac{C}{N}$ to estimate the volume of intersection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To approximate the volume of the intersection, we use:\n",
    "\n",
    "$$ V_N = V_D \\frac{1}{N} \\sum_{i=1}^{N} f(X_i) =  V_D \\frac{C}{N} $$\n",
    "\n",
    "where $ V_D$ is the volume of the domain used to generate the sample points. In this example, we considered the domain as the box around the intersection of the cylinders, hence \n",
    "\n",
    "$$ V_D = (2r)^3 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use your function `prob_inside_volume` to approximate the volume $V_{N}$ for $N = 1000$ for cylinders of radius 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 1000\n",
    "r = 1\n",
    "Vn = ...\n",
    "print(Vn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### d) Comparing with the exact solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two cylinders of radius r oriented long the z- and x-axes gives the equations $x^2+y^2=r^2$ and $y^2+z^2=r^2$ \n",
    "\n",
    "The volume common to two cylinders was known to Archimedes and the Chinese mathematician Tsu Ch'ung-Chih, and does not require calculus to derive. Using calculus provides a simple derivation, however. The volume is given by\n",
    "\n",
    "$$𝑉 = \\int_{-r}^{r}(2 \\sqrt{𝑟^2−𝑧^2})^2 𝑑𝑧= \\frac{16}{3}𝑟^3$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use your function `prob_inside_volume` to approximate the volume $V_{N}$ for increasing values of $N$ defined in `Nvalues`. Store each $V_{N}$ in a list `approxVol`. Plot $N$ vs $V_{N}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Nvalues = [(10**N) for N in range(1,7)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "approxVol = []\n",
    "\n",
    "# Add code here\n",
    "\n",
    "plt.plot(Nvalues,approxVol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the absolute error, using the exact expression given above. Plot $N$ vs $error$.  Compare with the known asymptotic behavior of the error $O(1/\\sqrt{N})$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "r = 1\n",
    "trueVol = (16.0/3.0)*r**3\n",
    "plt.loglog(Nvalues,np.abs(np.array(approxVol)-trueVol))\n",
    "plt.loglog(Nvalues, 1/np.sqrt(Nvalues), 'm')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random Walk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we will generate a list of random numbers and plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "series = np.random.rand(200)\n",
    "plt.plot(series)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A **random walk** is different from a list of random numbers because the next value in the sequence is a modification of the previous value in the sequence. Here is simple model of a random walk:\n",
    "\n",
    "- Start with a random number of either -1 or 1.\n",
    "- Randomly select a -1 or 1 and add it to the observation from the previous time step.\n",
    "- Repeat step 2 for as long as you like.\n",
    "\n",
    "Create the array `rand_walk` with N = 1000 points using the algorithm above.  Plot the array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 1000\n",
    "rand_walk = [0]\n",
    "\n",
    "# add code here\n",
    "\n",
    "plt.plot(rand_walk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 5) Modeling stock prices (simple model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that we are interested in investing in a specific stock. We may want to try and predict what the future price might be with some probability in order for us to determine whether or not we should invest.\n",
    "\n",
    "\n",
    "The simplest way to model the movement of the price of an asset in a market with no moving forces is to assume that the price changes with some random magnitude and direction, following a random walk model. \n",
    "\n",
    "$$p_{t} =  p_{t-1} + \\delta p $$\n",
    "\n",
    "We model the magnitude of the price change with the roll of a die. The function `dice` \"rolls\" an integer number from 1 to 6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dice():\n",
    "    return np.random.randint(1,7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To model prices increasing or decreasing, we will use the \"flip\" of a coin. The function `flip` returns either $-1$ (decreasing price) or $1$ (increasing price)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def flip():\n",
    "    return np.random.choice([1, -1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "By combining these two functions, we are able to obtain the price change at a given time. Here we will assume that a coin flip combined with a dice roll gives the price change for a given day."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### a) Performing one numerical experiment:\n",
    "\n",
    "Use the random walk model described above to predict the asset price for each day over a period of $N$ days. The initial price of the stock is $p0$. Store the price per day in the numpy array `price`.\n",
    "\n",
    "For now, use `N = 1000` and `p0 = 100`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 1000\n",
    "p0 = 100\n",
    "price = [p0]\n",
    "\n",
    "\n",
    "# Add code here\n",
    "\n",
    "\n",
    "plt.plot(price)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Does this plot resemble the short-term movement of the stock market?\n",
    "\n",
    "**Observations**:\n",
    "\n",
    "Performing one time step per day may not be enough to fully capture the randomness of the motion of the market. In practice, these N steps would really represent what the price might be in some shorter period of time (much less than a whole day).\n",
    "\n",
    "Furthermore, performing a single numerical experiment will not give us a realistic expectation of what the price of the stock might be after a certain amount of time since the stock market with no moving forces consists of random movements.\n",
    "\n",
    "Run the code snippet above several times (just do shift-enter again and again). What happens to the asset price after N days?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will be running several numerical experiments that simulates the price variation over N days. Wrap the code snippet above in the function `simulate_asset_price` that takes `N` and `p0` as argument and returns the array `price`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simulate_asset_price(N,p0):\n",
    "\n",
    "    \n",
    "    \n",
    "    \n",
    "    return np.array(price)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### b) Performing M=10 different numerical experiments, each one with N = 1000 days"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each numerical experiment, use the function `simulate_asset_price` with `N = 1000` and `p0 = 200`.\n",
    "\n",
    "Store all the M=10 arrays `price` in the 2d array `prices_M` with shape `(N,M)`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 1000 # days\n",
    "M = 10   # number of numerical experiments\n",
    "p0 = 200 # initial asset price"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "price_M = []\n",
    "for i in range(M):\n",
    "    price = simulate_asset_price(N,p0)\n",
    "    price_M.append(price)\n",
    "price_M = np.array(price_M).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then you can plot your results using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(price_M);\n",
    "plt.title ('M numerical experiments');\n",
    "plt.xlabel('Day');\n",
    "plt.ylabel('Price');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a more insightful prediction as to what the price of a given stock might be in the future. Suppose we want to predict the asset price at day 1000. We can just take the last element of the numpy array `price`!\n",
    "\n",
    "Create the variable `predicted_prices` to store the predicted asset prices for day 1000 for all the M=10 numerical experiments. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "predicted_prices = price_M[-1,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the histogram of the predicted price:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.hist(predicted_prices);\n",
    "plt.title('Asset price distribution at day 1000 from M numerical experiments')\n",
    "plt.xlabel('Asset prices')\n",
    "plt.ylabel('Number of Occurrences')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Go back and change the number of numerical experiments**. Set M = 1000 and run again. Better right?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can calculate the mean of the distribution to get the “expected value” for the stock on day 1000. What do you get?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "predicted_prices.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is one problem with our simple model. Our model does not incorporate any information about our specific stock other than the starting price. In order for us to get a more accurate model, we need to find a way incorporate the previous price of the stock. In the next example, we explore another model for stock price behavior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 6) Modeling stock prices (Black-Scholes Model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now model stock price behavior using the [Black-Scholes model](https://en.wikipedia.org/wiki/Black_Scholes_model), which employs a type of log-normal distribution to represent the growth of the stock price. Conceptually, this representation consists of two pieces:\n",
    "\n",
    "a) Growth based on interest rate only\n",
    "\n",
    "\n",
    "b) Volatility of the market\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Stock prices evolve over time, with a magnitude dependent on their volatility. The Black Scholes model treats this evolution in terms of a random walk. To use the Black-Scholes model we assume:\n",
    "\n",
    "   - some volatility of stock price.  Call this $\\sigma$\n",
    "   - a (risk-free) interest rate called $r$; and\n",
    "   - the price of the asset is [geometric Brownian motion](https://en.wikipedia.org/wiki/Geometric_Brownian_motion), or in other words the log of the random process is a normal distribution.\n",
    "\n",
    "which leads to the following expression for the predicted asset price:\n",
    "\n",
    "$$ S_T = S_0 e^{(r - \\frac{\\sigma^2}{2})T + \\sigma \\sqrt{T} \\,\\epsilon}$$\n",
    "\n",
    "meaning $ S_T/S_0 $ are normally distributed with mean $(r - \\frac{\\sigma^2}{2})T$ and variance  $\\sigma^2 T$, where \n",
    "\n",
    "   - $\\sigma$ is the volatility, or standard deviation on returns.\n",
    "   - $\\epsilon$ is a random value sampled from the normal distribution $\\mathcal{N}(0,1)$\n",
    "   - $S_T$ price of the asset at time $T$\n",
    "   - $S_0$ initial price of the asset \n",
    "   - $r$ is the interest rate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To predict the asset price at time $T$, we will discretize the total time to maturity in small steps $\\Delta t$. For each increment, we will use:\n",
    "\n",
    "$$ S_{t+\\Delta t} = S_t e^{(r - \\frac{\\sigma^2}{2})\\Delta t + \\sigma \\sqrt{\\Delta t} \\,\\epsilon}$$\n",
    "\n",
    "For example, if we want to obtain the asset price after 30 days, and we use the assumption that prices change with the increment of one day, then our total number of price estimates $S_t$ is $N = 30$, and $\\Delta t = 1/N$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Write the function `St_GBM` that will compute the price of an asset at time $t+\\Delta t$ given the parameters ($S_t,r,\\sigma,\\Delta t$)\n",
    "\n",
    "```python\n",
    "def St_GBM(St, r, sigma, deltat):\n",
    "    St_next = ... # Calculate this\n",
    "    return St_next\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This model now gives us a more accurate way to predict the future price.\n",
    "\n",
    "Assume that at the initial time $t = 0$ the asset price is $S0 = 100$, the interest rate is $r = 0.05$ and volatity is $\\sigma = 0.1$. \n",
    "\n",
    "Calculate the daily price movements using `St_GBM` for a period of 252 days (typical number of trading days in a year). Store your results in the array `price`. Then plot your results using `plt.plot(price)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 252 # number of days\n",
    "S0 = 100\n",
    "r = 0.05\n",
    "sigma = 0.1\n",
    "deltaT = 1/N\n",
    "\n",
    "\n",
    "# write code here\n",
    "\n",
    "\n",
    "plt.plot(price)   \n",
    "plt.title('Simulation of Price using Black-Scholes Model')\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Price')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Unfortunately volatility is usually not this small.... Run the code snippet above to predict the price movement for a volatity $\\sigma = 0.5$**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have managed to successfully simulate a year’s worth of future daily price data. Unfortunately this does not provide insight into risk and return characteristics of the stock as we only have one randomly generated path. The likelyhood of the actual price evolving exactly as described in the above charts is pretty much zero. We should modify the above code to run multiple numerical experiments (or simulations). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Perform M=10 different numerical experiments, each one with N = 252 days"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each numerical experiment, determine the array `price` using N = 252 days. Make sure to store all the M=10 arrays `price` in the 2d array `prices_M`.\n",
    "\n",
    "For this sequence of numerical experiments, assume that the initial asset price is `S0 = 100`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 252  # days\n",
    "M = 10000  # number of numerical experiments\n",
    "S0 = 100 # initial asset price\n",
    "r = 0.05\n",
    "sigma = 0.5\n",
    "deltaT = 1/N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then plot the result using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(price_M);\n",
    "plt.title ('M numerical experiments of Black-Scholes Model');\n",
    "plt.xlabel('Day');\n",
    "plt.ylabel('Price');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spread of final prices is quite large! Let's take a further look at this spread. Create the variable `predicted_prices` to store the predicted asset prices for day 252 (last day) for all the M=10 numerical experiments. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "predicted_prices = price_M[-1,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the histogram of the predicted prices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.hist(predicted_prices,30);\n",
    "plt.title('Predicted asset price distribution at day 252 from M numerical experiments')\n",
    "plt.xlabel('Asset prices')\n",
    "plt.ylabel('Number of Occurrences')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the mean and standard deviation of the distribution for the stock on the last day. What do you get?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Congratulations! You now have a prediction for a future price for a given stock."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 7) An example of a 2-d random walk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mariastenchia was in urgent need to use the restroom. Luckly, she saw Murphy's pub open and decided to go in for  relief. \n",
    "\n",
    "Unfortunately, Mariastenchia is not feeling very well, and due to some unknown reasons, she is confused and dizzy, and hence not sure if she can make it to the bathroom. After a quick evaluation, she decided that if she cannot get there in less than 300 steps, she will be in serious trouble.\n",
    "\n",
    "Do you think she can make a successful trip to the restroom? Let's help her estimating her odds.\n",
    "\n",
    "![](figures/murphy2.jpg)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " The helper function below plots the floor plan."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Helper function to draw the floor plan\n",
    "# You should not modify anything here\n",
    "def draw_murphy(wc,person,room_size):\n",
    "    fig = plt.figure(figsize=(6,6))\n",
    "    ax = fig.add_subplot(111)\n",
    "    plt.xlim(0,room_size)\n",
    "    plt.ylim(0,room_size)\n",
    "    plt.xlabel(\"X Position\")\n",
    "    plt.ylabel(\"Y Position\");\n",
    "    ax.set_aspect('equal')\n",
    "    \n",
    "    rect = plt.Rectangle(wc[:2], wc[-1], wc[-1], color=(0.6,0.2,0.1) )\n",
    "    ax.add_patch(rect)\n",
    "    plt.text(wc[0],wc[1]+wc[2]+0.2,\"WC\")\n",
    "\n",
    "    rect = plt.Rectangle((0,0),2,0.1, color=(0,0,1) )\n",
    "    ax.add_patch(rect)\n",
    "    plt.text(0.5,0.2,\"door\")\n",
    "   \n",
    "    circle = plt.Circle(person[:2], 0.3, color=(1,0,0))\n",
    "    ax.add_patch(circle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Let's take a look at the floor plan of Murphy's pub"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will simplify the floor plan with a square domain of size `room_size = 10`. The bottom left corner of the room will be used as the origin, i.e. the position `(x,y) = (0,0)`.\n",
    "\n",
    "The bathroom location is indicated by a square, with the left bottom corner located at `(8,8)` and dimension `h = 1`. These quantities are stored as a tuple `bathroom = (8,8,1)`.\n",
    "\n",
    " The street door is located at the bottom half, indicated by the blue line. Mariastenchia initial position is given by `initial_position = (1,1)`, marked with a red half-dot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "room_size = 10\n",
    "bathroom = (8,8,1) \n",
    "initial_position = (1,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will simplify Mariastenchia's challenge and remove all the tables and the bar area inside Murphy's. Here is how the room looks like in our simplified floor plan:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "draw_murphy(bathroom,initial_position,room_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How are we going to model Mariastenchia's walk?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Since Mariastenchia is dizzy and confused, we will model her steps as a random walk. \n",
    "\n",
    "- Each step will be modeled by a magnitude and direction in a 2d plane. We will assume the magnitude as 1.\n",
    "\n",
    "- The direction is given by a random angle $\\alpha$ between $0$ and $2\\pi$.\n",
    "\n",
    "- Combining the angle and magnitude, her step is defined as:\n",
    "\n",
    "$$ step = [\\cos(\\alpha), \\sin(\\alpha)]$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write the function `random_step` that takes as argument the current position, and returns the new position based on a random step with orientation $\\alpha$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def random_step(current_position):\n",
    "\n",
    "    # add code here\n",
    "    return(new_position)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make Mariastenchia give her 100 steps, using the function `random_step`. Complete the code snippet below, and plot her path from the door (given by the variable `initial_position` above) to her final location. Did she reach the bathroom?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 300\n",
    "position = [list(initial_position)]\n",
    "\n",
    "# update the list position\n",
    "    \n",
    "draw_murphy(bathroom,initial_position,room_size)\n",
    "x,y = zip(*position)\n",
    "plt.plot(x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You probably noticed Mariastenchia hitting walls, or even walking through them! Let's make sure we impose this constraints to the random walk, so we can get some feasible solution. Here are the rules for Mariastenchia's walk:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- If Mariastenchia runs into a wall, the simulation stops, and we call it a \"failure\" (ops!)\n",
    "- If the sequence of steps takes Mariastenchia to the restroom, the simulation stops, and we call it a success (yay!).\n",
    "- To simplify the problem, the \"restroom\" does not have a door, and Mariastenchia can \"enter\" the square through any of its sides.\n",
    "- Mariastenchia needs to reach the restroom in less than 100 steps, otherwise the simulation stops (at this point, it is a little too late for her...). This is also called a failure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `check_rules` checks if the `new_position` is a valid one, according to the rules above. The function returns `0` if the `new_position` is a valid step (still inside Murphy's and searching for the restroom), `1` if `new_position` is inside the restroom (sucess), and `-1` if `new_position` is a failure (Mariastenchia became a super woman and crossed a wall)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def check_rules(room_size,wc,current_position):\n",
    "    x,y,h = wc\n",
    "    # Checking if inside the room:\n",
    "    if ( (current_position[0] > 0) & (current_position[0] < room_size) & (current_position[1] > 0) & (current_position[1] < room_size)):\n",
    "        # Checking if found the restroom\n",
    "        if ( (current_position[0] > x) & (current_position[0] < x + h) & (current_position[1] > y) & (current_position[1] < y + h)):\n",
    "            return 1 \n",
    "        else:\n",
    "            return 0\n",
    "    else:\n",
    "        return (-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Modify the code snippet above, so that for every step, you check for the constraints using `check_rules`. Instead of giving all the 100 steps, you should stop earlier in case she reaches the restroom (`check_rules == 1`) or she hits a wall ( `check_rules == -1`) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 300\n",
    "position = [list(initial_position)]\n",
    "\n",
    "\n",
    "# add code here\n",
    "    \n",
    "draw_murphy(bathroom,initial_position,room_size)  \n",
    "x,y = zip(*position)\n",
    "plt.plot(x,y)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like this random walk does not give her much of a chance to get to the restroom. She may need more steps, or we can modify her walk to be a little less random. All you need to do is to modify your function `random_step`. What about we make her move forwards (in the positive direction of y) with a 70% probability (meaning she would move backwards with 30% probability). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def random_step_prob(current_position):\n",
    "\n",
    "    # add theta calculation\n",
    "    \n",
    "    new_position = current_position.copy()\n",
    "    new_position[0] += np.sin(theta*np.pi/180)\n",
    "    new_position[1] += np.cos(theta*np.pi/180)\n",
    "    return(new_position)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Let's estimate the probability Mariastenchia reaches the restroom:\n",
    "\n",
    "You should now run many simulations (one attempt to reach the restroom), and tally how many attempts are successful. The probability to reach the bathroom is `n_success/M`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "success = 0\n",
    "track_paths = []\n",
    "\n",
    "M = 100000\n",
    "N = 300\n",
    "\n",
    "for i in range(M):\n",
    "\n",
    "    position = [list(initial_position)]    \n",
    "\n",
    "    # One numerical experiment to check if successful trip to restroom or not\n",
    "            \n",
    "    if ((result == 1) & (not i%(M/100))):\n",
    "        track_paths.append(position)\n",
    "            \n",
    "\n",
    "draw_murphy(bathroom,initial_position,room_size)  \n",
    "for l in range(len(track_paths)):\n",
    "    x,y = zip(*track_paths[l])\n",
    "    plt.plot(x,y)\n",
    "\n",
    "print(\"probability is \", success/M)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "MonteCarloMethods.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
